C6.W2 NYT: Geodata processing¶
Performed by Alexander Manushin 10.02.2019
import pandas as pd
import numpy as np
from matplotlib import pylab as plt
from matplotlib import cm
import os
os.environ['PROJ_LIB'] = r'C:\Users\aaman\.conda\pkgs\proj4-5.2.0-hc56fc5f_1\Library\share'
from mpl_toolkits.basemap import Basemap
%matplotlib inline
1. Read data and sum trips¶
Загрузите агрегированные данные о поездках в мае 2016. Просуммируйте общее количество поездок такси из каждой географической зоны и посчитайте количество ячеек, из которых в мае не было совершено ни одной поездки.
Load aggregated data of trips for May 2016. Sum total trips' quantity from each geographical zone and count the number of cells of which not a single trip was made in May.
data = pd.read_csv('aggr_data.csv')
data.head()
data.groupby('pickup_region').sum()
print 'Zero trips regions: {}'.format(2500 - data.groupby('pickup_region').sum().shape[0])
2. Static map with ESB point¶
Нарисуйте статическую карту Нью-Йорка. Поставьте на карте точку там, где находится Эмпайр-Стейт-Билдинг.
Plot static map of New York City. Point the Empire State Building on map.
west, south, east, north = -74.25559, 40.49612, -73.70001, 40.91553
#creating map
m = Basemap(projection='merc', llcrnrlat=south, urcrnrlat=north,
llcrnrlon=west, urcrnrlon=east, lat_ts=south, resolution='f', epsg=2262)
fig = plt.figure(figsize=(10,10))
m.arcgisimage(service='ESRI_Imagery_World_2D', xpixels = 1500, verbose= True)
m.drawcoastlines()
plt.title('New York City map')
# Map (long, lat) to (x, y) for plotting
x, y = m(-73.985667, 40.748444)
plt.plot(x, y, 'ok', markersize=5, color='red')
plt.text(x, y, 'Empire State Building', fontsize=14, color='white')
3. Visualization on a map¶
Поверх статической карты Нью-Йорка визуализируйте данные о поездках из каждой ячейки так, чтобы цветовая шкала, в которую вы окрашиваете каждую ячейку, показывала суммарное количество поездок такси из неё.
Over the NYC static map visualize data of trips from each cell so that color scale that paints each cell shows the total trips' quantity from the cell.
#preprocessing of data
data['x']=(data.pickup_region.values)/50+1
data['y']=(data.pickup_region.values - 1)%50 + 1
trip_sum = pd.DataFrame(np.zeros([50,50]), index=range(1,51), columns=range(1,51))
pivot_sum = pd.pivot_table(data, values='trips', index=['y'], columns=['x'], aggfunc=np.sum).fillna(0)
trip_sum = trip_sum.add(pivot_sum)
trip_sum.replace(0, np.nan, inplace = True)
#matrix of trip sums
trip_sum
# map
fig = plt.figure(figsize=(10,10))
m.arcgisimage(service='ESRI_Imagery_World_2D', xpixels = 1500, verbose= True)
m.drawcoastlines()
plt.title('New York City map')
# ESB point
x, y = m(-73.985667, 40.748444)
plt.plot(x, y, 'ok', markersize=5, color='red')
plt.text(x, y, 'Empire State Building', fontsize=14, color='white')
# colorgrid
lons, lats = m.makegrid(50, 50)
x,y = m(lons,lats)
my_cmap = plt.get_cmap('Reds')
my_cmap.set_under('white')
cs = m.pcolormesh(x, y, trip_sum.fillna(0), cmap = my_cmap, alpha=0.4, zorder=2)
m.colorbar(cs, extend = 'min')
plt.show()
4. Interactive map with Statue of Liberty point¶
Вставьте интерактивную карту Нью-Йорка — такую, которую можно прокручивать и увеличивать. Поставьте метку там, где находится статуя свободы.
Add interactive NYC map - one that can be scrolled and zoomed. Point the Statue of Liberty.
# import folium interactive map
import folium
#plot map
m = folium.Map(location=[40.730610, -73.935242], tiles='OpenStreetMap', zoom_start=11)
folium.Marker([40.689247, -74.044502], tooltip='Statue of Liberty').add_to(m)
m
5. Interactive map with cells mean month trips¶
Нарисуйте на интерактивной карте Нью-Йорка ячейки так, чтобы их цвет показывал среднее за месяц количество поездок такси в час из этой зоны.
On the interactive NYC map plot cells so that their color shows average monthly quantity of taxi trips per hour from the zone.
#read regions data
regions = pd.read_csv('regions.csv', ';')
regions.head()
#preprocessing of geo data (creating poligons)
from geojson import Polygon, GeometryCollection, Feature, FeatureCollection, dump
def polygons(regions):
polygons = []
for i in range(len(regions)):
r_coord = [(regions.west.values[i],regions.south.values[i]),(regions.east.values[i],regions.south.values[i]),
(regions.east.values[i],regions.north.values[i]),(regions.west.values[i],regions.north.values[i]),
(regions.west.values[i],regions.south.values[i])]
f = Feature(geometry = Polygon([r_coord]), properties = {'reg_id':str(regions.region.values[i])})
polygons.append(f)
return FeatureCollection(polygons)
geo_data = polygons(regions)
#preprocessing of trips data
trips_mean = data.groupby('pickup_region').mean()
trips_mean.reset_index(inplace=True)
trips_mean.pickup_region = trips_mean.pickup_region.astype(str)
#plot a map
m = folium.Map(location=[40.730610, -73.935242], tiles='OpenStreetMap', zoom_start=11)
#plot SoL point
folium.Marker([40.689247, -74.044502], tooltip='Statue of Liberty').add_to(m)
#plot colorgrid
folium.Choropleth(
geo_data=geo_data,
name='choropleth',
data=trips_mean,
columns=['pickup_region', 'trips'],
key_on='feature.properties.reg_id',
fill_color='Reds',
nan_fill_color='gray',
fill_opacity=0.4,
line_opacity=0.1,
legend_name='Mean trips (%)'
).add_to(m)
m
6. Checking the regions from which the trip is impossible¶
Чтобы не выбирать из всех 2500 ячеек вручную, отфильтруйте ячейки, из которых в мае совершается в среднем меньше 5 поездок в час. Посчитайте количество оставшихся. Проверьте на карте, что среди этих ячеек нет таких, из которых поездки на самом деле невозможны.
In order not to select manually from all 2500 cells, filter out cells from which on average less than 5 trips per hour are made in May. Count the number remaining. Check on the map that among these cells there are none of which trips are actually impossible.
#data filtration
trips_mean_filtered = trips_mean[trips_mean.trips > 5]
trips_mean_filtered.sort_values('trips').head()
trips_mean_filtered.info()
#save trips into csv
trips_mean_filtered.to_csv('trip_mean_filtered.csv', index=False)
print '{} regions has more then 5 trips per hour'.format(len(trips_mean_filtered))
#plot a map with filtered data
m = folium.Map(location=[40.730610, -73.935242], tiles='OpenStreetMap', zoom_start=11)
#plot SoL poin
folium.Marker([40.689247, -74.044502], tooltip='Statue of Liberty').add_to(m)
#plot grid
folium.Choropleth(
geo_data=geo_data,
name='choropleth',
data=trips_mean_filtered,
columns=['pickup_region', 'trips'],
key_on='feature.properties.reg_id',
fill_color='Reds',
nan_fill_color='gray',
fill_opacity=0.4,
line_opacity=0.1,
legend_name='Mean trips (%)'
).add_to(m)
m
Checking result: trips are possible from all regions which has more then 5 trips per hour
7. Jupyter notebook saving¶
Сохраните ноутбук в формате html, запакуйте ipynb и html в один архив и загрузите его в форму (html облегчит задачу проверки вашим рецензентам — в ipynb динамические карты часто не сохраняются, а для выполнения кода рецензентам придётся устанавливать все библиотеки).
Save the jupyter notebook in html format, pack ipynb and html into one archive and upload it to the form (html will make it easier for your reviewers to check - dynamic cards are often not saved in ipynb, and reviewers will have to install all libraries to execute the code).